00(o  -  S34) 


AD-A236  067 


*  4'  ^ 


l 


i 

v  > 


^  t  approved 

for  PUBLIC  RLJ.LASI 


91-01062 


A  One-Dimensional  Flux- Corrected 
Transport  Code  for 
Detonation  Calculations 

*  t 

D.A.  Jones,  E.S.  Oran  and  R.  Guirguis 


MRL  Research  Report 
MRL-RR-2-90 


Abstract 


loa  for 

DTIC  T*» 
'JotiaocvriortcJ 

Juatiricatlcxi. 


»T  — _ 

DXstrif  • 

Avaii-c.i ;.:r , 

,Avi\  1  i  ; 
I  ;  pen  1*1 


l«\-\ 


The  devel'jpment  of  a  one-dimensional  Flux- Corrected  Trans  (ml  code  to  model  detonation  in 
a  homogeneous  medium  is  described.  The  material  flow  is  modelled  using  the  Euler 
equations,  and  the  chemical  kinetics  by  a  two-step  induction  parameter  model  which  uses  a 
quasi-steady  induction  time  and  first  order  Arrhenius  kinetics.  Two  different  modes  of 
initiation  are  compared.  Conditions  necessary  for  a  self-sustaining  detonation  are  described 
and  illustrated.  A  detailed  comparison  is  made  between  the  variable  profiles  calculated  by 
the  code  and  those  calculated  aruilytically  using  the  simple  Chapman-Jouguet  theory  and 
self-similar  analysis,  and  the  overall  agreement  is  excellent.  The  effect  of  the  computational 
cell  size  on  these  solutions  is  also  considered. 
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A  One -Dimensional  Flux-Corrected 
Transport  Code  for  Detonation  Calculations 


1.  Introduction 

A  detonation  wave  is  a  supersonic  shuck  wave  travelling  through  a  reactive 
medium.  The  high  temperature  generated  by  the  passage  of  the  shock  through  the 
material  initiates  chemical  reaction  and  energy  release.  A  stable  detonation  with  a 
unique  detonation  velocity  will  develop  if  the  chemical  reaction  zone  can  remain 
coupled  to  the  wave  and  continuously  feed  energy  to  the  shock  front.  Numerical 
simulations  of  propagating  detonations  therefore  require  simultaneous  solutions  of 
equations  which  describe  both  the  material  flow  and  reaction  kinetics. 

Good  reviews  of  the  early  work  on  the  numerical  simulation  of  detonation  can  be 
found  in  the  books  by  Mader  [1]  and  Fickett  and  Davis  [2).  One  of  the  first 
calculations  was  reported  by  Hubbard  and  Johnson  [3].  Most  of  this  early  work 
relied  on  the  artificial  viscosity  method  of  von  Neumann  and  Richtmyer  [4]  to 
model  shock  formation,  although  one  exception  was  the  work  of  Fickett  and  Wood 
[5],  which  used  a  method-of-charactcristics  code  to  model  one-dimensional 
longitudinal  instabilities  in  overdriven  detonations.  Codes  which  are  based  on  the 
method-of-characteristics  technique  have  the  advantage  of  very  accurately  resolved 
shock  fronts,  but  they  are  not  very  easily  extended  into  two  dimensions. 

Conversely,  codes  based  on  the  artificial  viscosity  method  arc  readily  extended  to 
two  or  three  dimensions,  but  they  have  the  disadvantage  of  smearing  the  shock  front 
over  several  cells.  Artificial  viscosity  codes  also  cause  problems  when  coupled  to 
modem  schemes  which  model  the  reaction  kinetics  of  heterogeneous  condensed 
phase  explosives.  Many  of  these  schemes  incorporate  a  material  viscosity  term  to 
model  hot  spot  formation  [6],  and  the  identification  of  the  correct  value  of  this  term 
is  often  complicated  by  the  smearing  of  the  shock  front  caused  by  the  artificial 
viscosity.  This  problem  has  been  noted  in  a  recent  MRL  report  [7], 

Advances  in  computational  fluid  dynamics  during  the  1970s  have  removed  the 
need  for  reliance  on  the  artificial  viscosity  method  lor  two  or  three-dimensional 
shock  calculations.  The  Flux-Corrcctcd  Traasport  (FCT)  technique  developed  by 
Boris  and  Book  (8-  10J  showed  that  nonlinear  monotone  methods  could  accurately 


resolve  shock  fronts  to  within  one  cell  width  without  the  need  for  an  additional 
artificial  viscosity  term.  Since  then  a  variety  of  non-linear  methods  have  been 
derived  and  applied  extensively  to  gas  dynamic  and  fluid  calculations.  The  paper 
by  Sod  [11]  provides  a  good  review  of  some  of  these  methods.  In  general,  however, 
the  explosives  community  has  been  slow  in  applying  these  techniques  to 
the  numerical  simulation  of  detonation  phenomena.  One  exception  to  this  is  the 
reactive  flow  work  of  Oran  and  collaborators  at  the  Naval  Research  Laboratory. 
They  have  used  one  and  two-dimensional  FCT  codes  to  model  both  gaseous 
detonation  [12,  13]  and  homogeneous  condensed  phase  detonations  [14-16].  A 
comprehensive  review  of  this  work  has  recently  been  published  by  Oran  and  Boris 
[17],  Nunziato  and  Baer  have  also  used  FCT  methods  to  model  flame  propagation 
and  growth  to  detonation  in  multiphase  flows  [18],  and  Thomas  has  also  used  FCT 
codes  for  detonation  calculations  [19]. 

MRL  is  interested  in  the  shock  initiation  of  heterogeneous  condensed  phase 
explosives.  A  previous  MRL  report  [7]  has  critically  reviewed  models  for  the  shock 
initiation  of  these  materials  and  their  suitability  for  implementation  in  various 
hvdrocodes.  This  report  describes  a  one-dimensional  FCT  computer  code  which 
models  the  propagation  of  detonation  in  homogeneous  materials  with  a  polytropic 
equation  of  state.  The  decomposition  of  the  material  is  described  by  first-order 
kinetics  with  an  Arrhenius  temperature  dependence.  A  pressure  and  temperature 
dependent  induction  time  is  also  included.  Replacement  of  the  polytropic  equation 
of  state  with  one  more  suited  to  condensed  phase  materials,  and  implementation  of 
one  of  the  schemes  described  in  [7],  would  then  allow  the  code  to  model  shock 
initiation  of  heterogeneous  explosives.  The  code  is  similar  to  those  described  in 
broad  outline  in  [17],  and  in  more  detail  in  [14],  In  this  report  we  present  a  more 
detailed  description  of  some  aspects  of  the  model  and  its  numerical  solution.  We 
make  a  detailed  comparison  between  the  variable  profiles  calculated  by  the  code 
and  those  calculated  analytically  using  the  simple  Chapman-Jouguct  (CJ)  theory', 
and  we  also  consider  the  effect  of  the  computational  cell  size  on  these  solutions. 
Finally,  we  illustrate  the  use  of  the  code  by  discussing  its  application  to  recent  work 
on  detonation  in  a  stoichiometric  mixture  of  hydrogen  and  oxygen. 


2.  The  Physical  Model 

The  numerical  simulation  of  detonation  requires  the  simultaneous  solution  of  the 
coupled  equations  describing  material  flow  and  chemical  reaction.  The  pressures 
generated  by  the  detonation  of  a  condensed  phase  explosive  arc  so  high  (typically  a 
tew  hundred  kilobars)  that  material  strength  may  be  neglected,  and  the  appropriate 
equations  describing  the  material  flow  are  those  of  reactive  fluid  dynamics  [20]. 

We  also  make  the  usual  assumption  that  energy  transport  by  heat  conduction, 
viscosity,  and  radiation  is  negligibly  small  compared  with  transport  by  motion  [2], 
In  this  case,  the  appropriate  fluid  equations  are  the  Euler  equations,  a  set  of  three 
coupled  partial  differential  equations  describing  the  conservation  of  mass, 
momentum,  and  energy.  In  a  one -dimensional  cartesian  coo;du  k  system,  these 
arc 
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where  t  denotes  time,  x  the  spatial  coordinate,  p  the  density,  v  the  fluid  (or  particle) 
velocity,  and  P  the  pressure.  The  quantity  E  is  the  total  energy  per  unit  volume 
and  is  defined  by 


E  -  p  e  *  0.5  p  v2 


(4) 


where  e  is  the  specific  internal  energy  (including  the  stored  chemical  energy). 

The  thermodynamic  properties  of  a  nonreacting  material  are  completely 
specified  by  the  fundamental  equation,  which  expresses  the  internal  energy  of  the 
system  in  terms  of  its  entropy,  volume,  and  mole  numbers.  The  partial  derivatives 
of  the  fundamental  equation  with  respect  to  its  dependent  variables  then  define  the 
equations  of  state  of  the  material  [21  ].  The  thermodynamic  state  can  also  be 
specified  precisely  if  all  the  equations  of  state  of  the  system  arc  known.  For  a 
simple  single-component  system  there  are  three  equations  of  state.  Only  two  of 
these  need  to  be  specified  however  as  the  Gibbs-Duhem  relation  may  always  be 
used  to  obtain  the  third  equation  [2 1 J.  The  two  equations  of  slate  normally 
specified  are  the  mechanical  equation  of  state  and  the  caloric  equation  of  state, 
which  have  the  forms 


P  =  P(e,  p) 


T  =  T(e,  p) 


where  T  is  the  temperature  of  the  material.  We  use  a  polytropic  model  for  our 
fluid,  which  describes  an  ideal  gas  for  which  the  ratio  of  specific  heats  at  constant 
pressure  and  temperature  is  a  constant.  The  equations  of  state  are  then  given  by 


[221, 


P  -  ( y  -  l)pe 


(5) 


T  - 


e 


(6) 


where  y  =  c  ic  ,  and  c  and  c  arc  the  specific  heals  at  constant  pressure  and 
*  p  v  p  v  . 

volume  respectively.  Equations  (5)  and  (6)  can  also  be  combined  to  give 


P  -  pRT 


(7) 


in  which  R  is  the  specific  gas  constant,  R  =  R  irriw.  where  R  is  the  universal  gas 
constant  and  m\v  is  the  molecular  weight  of  the  fluid. 

The  material  is  allowed  to  undergo  a  single  irreversible  chemical  reaction, 
represented  schematically  by 


explosive  reactants  — »  explosive  products 


(8) 


The  condition  that  both  die  molecular  weight  and  '/remain  unaltered  during  die 
progress  of  the  reaction  is  also  imposed.  The  chemical  composition  at  any  time  is 
described  by  the  progress  variable  u  .  which  we  define  to  be  the  explosive  mass 
fraction,  i.e. 

mass  of  reactants  (9) 

mass  of  reactants  +  mass  of  products 


w  therefore  goes  from  one  to  zero  as  the  reaction  proceeds  from  pure  reactants  to 
pure  products.  The  temperature  increase  as  the  reaction  proceeds  is  calculated  from 


T  (c  -  wAE) 


(10) 


where  A E  is  the  heat  of  reaction.  The  pressure  is  then  calculated  from  equation  (7). 
This  is  a  commonly  used  energetic  scheme  in  detonation  calculations  |23|. 

The  chemical  reaction  represented  schematically  by  equation  (8)  often  takes 
place  in  two  distinct  stages.  In  the  first  step,  the  explosive  fuel  is  converted  into 
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radicals  by  a  series  of  endothermic  processes.  During  this  step,  the  concentration 
of  radicals  steadily  increases  but  remains  so  small  that  the  amount  of  fuel  consumed 
and  products  produced  can  be  neglected,  and  there  is  no  apparent  rise  in  either 
temperature  or  pressure.  We  allow  for  this  induction  period  by  the  explicit 
introduction  of  an  induction  time  x°  Induction  times  arc  measured  at  constant 
temperature  and  pressure,  and  we  use  a  superscript  zero  to  denote  that  x°  is  the 
induction  time  measured  (or  calculated)  under  these  conditions.  In  the  course  of  our 
calculations,  however,  a  given  control  volume  may  undergo  a  scries  of  temperature 
and  pressure  changes  and  the  true  induction  time  x  will  not  necessarily  be  equal  to 
x°  Therefore  x  is  determined  from  the  solution  of  the  integral  equation 


r  a!  -  1  (ID 

Jo  x  °(T,P) 

This  formula  provides  a  convenient  way  of  averaging  over  changing  temperature 
and  pressure  during  the  chemical  induction  period,  and  reduces  automatically  to  the 
correct  induction  delay  when  the  pressure  and  temperature  are  constant. 

We  also  introduce  an  induction  parameter  /.  which  represents  the  fraction  of 
induction  time  al reads'  elapsed.  >  is  obtained  from  the  solution  of  the  equation 


v. _ L_ 

dl  x°(T,P) 


(12) 


where  fun  =  <  Hnergy  release  is  then  initiated  when  f  -  I  A  useful  expression 
lor :  which  has  been  used  previous!}  |  lh)  has  the  lomi 


\°  (T,P) 


exp  (EZ/RT) 


r 


(13) 


where  P  is  a  reference  pressure. 

When  the  induction  time  has  elapsed,  the  energy  release  is  assumed  to  follow 
first  order  kinetics  with  an  Arrhenius  temperature  dependence 

tL  -  -  *Arcxp(-Er/RT),  d4> 

di 

where  Ef  is  the  activation  energy  and  A  is  the  frequency  factor.  Values  for  the 
constants  A  E  and  Af.  Er  can  be  found  from  solutions  of  the  detailed  chemical 
kinetic  equations  if  the  full  reaction  scheme  is  known  ( 13],  or  else  fitted  to 
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experimental  measurements  [15].  The  combination  ol  an  induction  parameter 
model  of  the  form  of  equation  (12)  and  reaction  rale  dependence  given  by  equation 
(14)  was  first  described  by  Oran  et  al.  [  24  ] . 

For  some  reactive  compositions,  both  the  induction  time  and  reaction  rate  arc 
independent  of  temperature  and  pressure  in  the  region  of  interest  and  can  be  taken 
as  constants.  We  describe  one  example  of  this  in  a  later  section  where  we  develop 
a  simple  model  for  the  detonation  of  a  stoichiometric  mixture  of  hydrogen  and 
oxygen. 

When  the  reaction  occurs  in  a  fluid  moving  with  a  velocity  v,  the  time 
derivatives  in  equations  (12)  and  ( 14)  denote  a  derivative  following  the  fluid  How 
and  the  equations  become 


8/+v8/.  1 

Tt  "51  7 


(15) 


and 


r)  H'  i; 

H T  +  v17 


w/trexp(  -  Er/RT) 


(16) 


respectively.  F.quations  (15)  and  1 16)  can  then  be  combined  with  the  continuity 
equation,  equation  ( 1  >.  to  give 


3(p f)  ,  3(p/v)  _  _P_  (17) 

ft  x " 


and 


8<P>y).  *  a(Pw)  -  -  p  nv4,.exp(  -  Er/RT)  tf«) 

dt  dx 


Equations  ( 17)  and  (18)  are  now  in  conservation  form  and  can  be  solved  by  the 
same  method  used  to  solve  equations  Cl)  through  (3). 

In  the  numerical  solution  procedure  outlined  in  the  next  section,  equation  ( 1 7)  is 
modified  slightly  to  convcct  the  product  w/  rather  than  /  alone,  i.c.  equation  (17) 
is  replaced  by 
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(19) 


3(p  wf)  d(pwfv)  wp 

—jr~ +  r* 


which  is  equivalent  to  writing  equation  (12)  in  the  form 


d(*f)  _  (20) 

dt  T° 


1  he  replacement  of  equation  (17)  by  equation  (19)  facilitates  calculations  in  cells 
containing  fractions  of  both  products  and  reactants.  Equation  (14)  ensures  that  any 
products  in  a  mixed  cell  have  no  effect  on  the  reaction  of  the  remaining  unreacted 
explosive.  The  reasoning  behind  the  replacement  cf  equation  (17)  by  equation  (19) 
is  complex,  but  is  explained  in  detail  in  reference  [15]. 


3.  Numerical  Solution 

The  computer  code  RSHOCK  solves  equations  (1),  (2),  (3),  (18)  and  (19)  using  the 
FCT  subroutine  JPBFCT  (a  computer  listing  of  RSHOCK  may  be  obtained  from 
D  A.  Jones  on  request).  Operator  splitting  is  used  for  the  chemical  source  terms  in 
equations  ( 1 8)  and  1 19),  and  equations  (4),  (7)  and  (10)  define  the  additional 
thermodynamic  v  ariables  in  terms  of  the  convccted  quantities.  JPBFCT  is  a 
slightly  modified  version  of  the  subroutine  ETBFCT,  which  is  described  in  detail  in 
reference  |25j.  This  multiple  entry'  subroutine  solves  generalised  continuity 
equations  of  the  form 


— .  4-  pv> 

r«  '  or 
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dD2 
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dp 

TT 


(21) 


using  an  algorithm  which  gives  fourth  order  accurate  phases  and  minimal  residual 
diffusion.  The  routine  can  treat  cartesian,  cylindrical  or  spherical  one-dimensional 
systems  for  a  =  1, 2,  3  respectively. 

The  solution  procedure  for  the  coupled  set  of  equations  is  to  first  call  JPBFCT  te 
convect  the  density  p.  Next,  the  subroutine  SOURCD  is  called  to  calculate  the 
source  term  for  the  momentum  equation,  and  then  JPBFCT  is  called  to  convect  the 
product  pv.  SOURCD  is  then  called  again  to  calculate  the  source  term  for  the 
energy  equation,  and  then  JPBFCT  called  to  advance  the  total  energy  £.  This 
procedure  is  performed  twice  for  each  hydrodynamic  lime  step  At.  On  th<1  first 
pass,  the  variables  p,  pv  and  £  are  advanced  over  a  time  0.5  Af,  and  are  then  used  to 
calculate  the  source  terms  at  this  intermediate  time.  On  the  second  pass  the 
variables  are  convected  over  the  full  step  At  using  the  source  term  values  evaluated 
at  the  half  step.  This  procedure  greatly  increases  the  accuracy  of  the  simulation  by 
correcting  a  phase  lag  in  the  values  of  the  variables  which  occurs  as  the  equations 
arc  solved  sequentially.  The  variables  pw  and  pwf  arc  also  convected  during  this 
sequence,  with  the  source  terms  on  the  right  side  of  equations  (18)  and  (19)  set  to 
zero.  Only  p  w  is  convected  on  the  half  step  cycle  because  the  pressure  depends  on 
w,  but  not/.  On  the  full  step  calculation,  both  vv  and/are  convected.  At  the 
completion  of  the  transport  stage,  e,  T  and  P  arc  then  calculated  from  equations  (4), 
(6)  and  (7). 

Next,  the  program  checks  to  see  if  the  induction  time  has  elapsed  by  calculating 
the  variable  DT,  which  is  equal  to  the  amount  of  time  remaining  before  /  equals 
one.  At  this  stage,  there  are  two  possibilities,  DT  is  either  greater  or  less  than  the 
hydrodynamic  time  step  Ar.  The  former  case  implies  that  no  burning  occurs  at  the 
end  of  the  lime  step,  and  so  the  time  dependent  part  of  equation  (19)  is  solved  via 
the  explicit  equation 


(H f)ne"  -  (H f)°ld  -At  WM/X 


(22) 


to  update  /,  and  w  remains  fixed  at  one.  The  latter  case  implies  that  burning  starts 
before  the  end  of  the  current  time  step.  The  quantity  DT  is  then  reset  to  the  bum 
time  within  the  time  step,  i.e.  DT  =  At  -  DT,  and  w  is  updated  by  the  implicit 
solution  of  the  time  dependent  part  of  equation  (18), 


new  _  w°U/(\ 


A,cxp(-E,/RT)DT) 


(23) 


Equations  (10)  and  (7)  are  then  used  to  update  the  temperature  and  pressure  after 
burning. 

Because  both  pressure  and  temperature  are  variables  which  are  not  convected  in 
conservation  form,  but  rather  are  derived  from  such  variables,  they  are  particularly 
prone  to  numerical  undershoots  and  overshoots.  Consequently,  at  various  places  in 
the  code  both  pressure  and  temperature  are  limited  so  that  unrealistic  undershoots, 
in  particular,  do  not  occur.  Similarly,  both  /  and  w  are  variables  which  are  only 
defined  between  the  limits  zero  and  one,  and  so  both  variables  are  limited  within  the 
code  to  ensure  that  these  conditions  apply.  In  particular,  conditions  are  imposed  so 
that  w  remains  fixed  at  one  while  /  is  less  than  one,  and  /  is  fixed  at  one  when  w 
is  less  than  one. 

The  hydrodynamic  time  step  At  is  determined  at  the  start  of  the  integration  step 
via  the  Courant  condition  which  we  use  in  the  form 


A i  -  0.25  min  ( 5 jc / ( I v I  +  c)] 


(24) 


where  8.v  is  the  computational  cell  size  and  c  is  the  sound  speed,  which  is  calculated 
from 


c:  -  yP/p 


(25) 


4.  Results  for  HJOJAr 

Our  first  application  of  the  code  is  to  a  stoichiometric  mixture  of  hydrogen  and 
oxygen  diluted  with  argon.  The  composition  of  the  system  is  H^Oj/Ar  in  the  ratio 
2:1:7.  This  system  has  been  studied  extensively  at  the  Naval  Research  Laboratory 
using  the  CHEMOD  code  and  a  detailed  description  of  the  full  kinetic  scheme  which 


converts  reactants  into  products  [26].  The  time  dependence  of  the  reaction  has 
been  followed  for  many  values  of  the  initial  temperature  and  pressure  and  the 
induction  time  and  rate  of  energy  release  have  been  found  to  closely  follow  the 
expressions  given  by  equations  (13)  and  (14),  respectively.  The  values  of  the 
constants  A^,  Ex,  Af  and  Er  found  by  this  procedure  are  given  in  Table  1.  Also 
shown  in  Table  1  are  the  vdues  of  the  remaining  constants  which  need  to  be 
specified  for  the  model,  i.e.  mw ,  y  and  A E.  As  we  are  using  a  model  in  which 
both  mw  and  y  remain  fixed  during  the  progress  of  the  reaction,  the  values  we  use 
for  these  constants  when  applied  to  a  reaction  in  which  mw  and  y  are  different  for 
reactants  and  products  need  some  consideration.  Our  only  constraint  should  be  that 
the  values  we  use  should  lie  somewhere  between  those  for  the  reactants  and  the 
products.  For  our  model  of  H^Oj/Ar,  we  have  chosen  to  calculate  mw  and  y  for 
the  reactant  composition.  Taking  the  molecular  weights  of  H2,  02  and  Ar  to  be  2, 
32  and  40,  respectively,  leads  to  a  molecular  weight  of  31 .6  for  the  reactant  mixture. 
The  y  for  the  reactants  is  calculated  from  the  expression 
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(26) 


and  leads  to  the  value  1.5555.  A E  is  defined  as  the  energy  release  per  unit  mass  of 
mixture.  The  formation  of  18  g  of  H20  releases  57.8  kcal  or  242  kJ  of  energy. 
Dilution  with  argon  though  results  in  an  energy  release  of  1 .53  kJ/g.  This  is  not  the 
value  used  for  A £  however  as  a  large  fraction  of  this  energy  will  not  be  available  to 
drive  the  detonation,  but  will  instead  lead  to  some  degree  of  dissociation  and 
ionisation  of  the  products.  We  assume  that  approximately  half  the  energy  released 
is  lost  to  these  processes,  and  so  wc  take  for  A E  a  value  of  0.75  kJ/g  (or 
0.178  keal/g). 

The  manner  in  which  a  steady  detonation  is  initiated  in  the  code  is  of  some 
importance.  In  the  results  we  present  for  HVOj/Ar  the  code  was  initiated  by 
depositing  an  excess  amount  of  energy  into  a  few  cells  adjacent  to  the  left  boundary 
of  the  grid.  The  exact  amount  of  energy,  and  the  number  of  cells  into  which  it 
should  be  deposited,  can  only  be  determined  by  trial  and  error.  A  handy  rule  of 
thumb  is  that  the  excess  energy  should  be  approximately  five  times  the  energy  in  the 
detonation  front.  For  a  computational  cell  size  of  4  x  10'4  m,  about  10  cells  were 
found  to  be  adequate  to  initiate  detonation.  In  all  the  calculations  presented  here, 
we  imposed  rigid  wall  boundary  conditions  at  either  end  of  the  grid. 
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Table  1: 

Parameter  Values  for 

:Ar  in  ratio  2:1 :7 

A 

r 

= 

24.0  x  108  s  1 

mw 

=  31.6 

Er 

= 

3.0  x  104  kcal 
(12.6  x  104  kJ) 

y 

-  1.5555 

At 

= 

5.684  x  10  8  s  1 

A E 

=  3.14  kcal/g'1 

(13.14  Ug- 

Ei 

= 

1.503  x  104  kcal 
(6.288  x  104  kJ) 

Figure  1  shows  simulation  results  for  this  model  after  20  000  time  steps.  The 
computational  grid  contained  4  500  cells  with  Ax  -  4.0  x  10  4  m  and  the  time  step 
had  an  average  value  of  0.05  ps.  The  detonation  was  initiated  by  depositing  an 
excess  energy  of  3.0  x  106  J  per  cell  into  the  first  10  cells.  The  initial  temperature 
and  pressure  were  300  K  and  1 .0  x  105  Pa  respectively,  and  the  run  took 
approximately  9  hours  of  VAX  8700  CPU  time.  The  pressure  profile  shows  a  very 
sharp  shock  front  and  a  jump  to  a  von  Neumann  point  of  about  21.5  x  105  Pa.  The 
von  Neumann  temperature  is  approximately  1840  K  and  estimates  of  the  induction 
time  and  rate  of  reaction  from  equations  (13)  and  (14)  using  these  values  indicate 
that  the  total  reaction  zone  should  span  about  six  computational  cells.  On  the  scale 
of  Figure  1 ,  this  means  that  the  fall  in  pressure  between  the  von  Neumann  point  and 
the  CJ  point  will  be  compressed  into  the  width  of  approximately  one  line,  and  this 
explains  the  sharp  spikes  displayed  by  the  pressure,  density  and  panicle  velocity 
plots,  as  well  as  die  abrupt  transition  between  zero  and  one  displayed  by  the  reaction 
variable  profile.  We  estimate  the  CJ  pressure  to  be  12.4  x  105  Pa,  and  following 
this,  we  observe  a  Taylor  wave  expansion  down  to  a  steady-state  pressure  of  4.5  x 
10- Pa. 

One  of  the  advantages  of  the  model  wc  have  described  is  that  analytical  results 
can  be  compared  to  the  simulation  results.  Using  simple  expressions  derived  in 
Thompson  135 1,  or  Landau  and  Lifshitz  [36],  we  calculate  a  CJ  velocity  for 
FT/Oj/Ar  of  1542  m/s.  Figure  1  also  shows  detonation  velocity  versus  time.  After 
5000  time  steps  (262  ps),  the  velocity  is  approximately  1570  m/s,  or  widiin  2%  of 
the  CJ  value.  After  a  further  15  000  time  steps  (at  l  =  1050  ps),  the  velocity  has 
decreased  to  1547  m/s,  which  is  0.3%  higher  than  the  CJ  value.  The  estimated  CJ 
pressure  of  12.4  x  105  Pa  is  1.6%  higher  than  the  analytical  value  of 
12.2  x  105  Pa,  while  the  estimated  CJ  temperature  of  2300  K  is  in  very  good 
agreement  with  the  analytical  value  of  2298  K.  A  fairly  large  (approximately  20%) 
temperature  overshoot  is  evident  in  the  plot  of  temperature  versus  computational 
cell  number.  Ideally  we  would  expect  the  temperature  to  be  shocked  to  the  von 
Neumann  point  at  1840  K,  and  then  rise  again  as  the  reaction  proceeds  to  a  peak 
value  of  2300  K  at  the  CJ  point.  Although  both  overshoots  and  undershoots  are  not 
unexpected  in  quantities,  the  overshoot  in  temperature  observed  here  is 
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Figure  I:  Simulation  results  for  HJOjAr  after  20  000  time  steps.  NDUMP  =  10, 
DSTRBO  =  3.0  x  106  J  and  A x  =  4.0  x  104  m.  The  filled  circles  on  the  pressure 
and  particle  velocity  profiles  are  from  the  self-similar  analysis  ofKuhl  and  Seizew 
[28j. 
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Figure  l  (continued):  Simulation  results  for  H^O^Ar  after  20  000  time  steps. 
NDUMP  =  10,  DSTRBO  =  3.0  x  IQ6  J  and  Ax  =  4.0  x  I  O'4  m.  The  filled  circles  on 
the  pressure  and  particle  velocity  profiles  are  from  the  self-similar  analysis  ofKuhl 
and  Seizew  [28/. 
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caused  by  the  narrowness  of  the  reaction  zone.  We  will  present  results  below 
which  she  *v  that  increasing  the  number  of  computational  cells  within  the  reaction 
zone  reduce?  this  overshoot.  The  peak  pressure  of  21.5  x  105  Pa  also  undershoots 
the  expectef  value  of  23.9  x  105  Pa  by  10%  but  this  is  not  related  to  the  resolution 
of  the  reaction  zone;  it  reflects  the  limitations  of  the  numerical  scheme  employed. 
The  other  feature  worthy  of  comment  is  the  very  high  temperature  at  the  left 
boundary  of  the  grid  which  is  still  present  after  20  000  time  steps.  This  is  a  remnant 
of  the  extra  energy  used  to  initiate  the  detonation  and  its  slow  decay  is  due  to  the 
absence  of  a  heat  diffusion  term  in  the  model,  as  well  as  the  zero  fluid  velocity 
boundary  condition.  We  will  shortly  present  a  different  method  of  initiation  which 
overcomes  this  problem. 

The  resolution  of  Figure  1  is  such  that  the  detonation  may  be  viewed  as  an 
idealised  CJ  detonation.  In  this  case  the  flow  field  is  self-similar  and  may  be 
calculated  using  established  techniques  [27],  Kuhl  and  Seizew  [28]  have  analysed 
ideal,  strong,  CJ  detonations  using  the  phase-plane  technique  and  presented 
tabulated  solutions  for  the  flow  variables  for  planar,  cylindrical  and  spherical 
geometries  for  detonations  in  gases  (y  =  1.3)  and  solids  (y  =  2.7).  All  properties  of 
the  flow  field  behind  strong  CJ  detonations  are  functions  of  the  detonation  velocity 
and  y  alone  (and  boundary  conditions).  Kuhl  and  Seizew  explored  the  effects  of 
geometry  and  changes  in  y  on  the  computed  flow  fields  and  found  that  pressure  and 
velocity  profiles  showed  only  a  weak  dependence  on  y,  while  density  and  internal 
energy  profiles  showed  a  pronounced  dependence.  Hence  we  can  use  their 
tabulated  solutions  for  gases  (y  =  1.30)  to  analyse  our  computed  pressure  and 
velocity  profiles  for  Hj/Oj/Ar  (y  =  1.5555).  The  results  are  shown  on  the  pressure 
and  velocity  profiles  in  Figure  1.  The  good  agreement  shows  that  the  1-D  code  has 
accurately  calculated  the  flow  field  behind  the  detonation  front. 

It  is  important  to  ensure  that  the  results  presented  above  are  independent  of  the 
computational  cell  size.  From  the  excellent  agreement  shown  between  the 
numerical  simulation  and  the  analytically  calculated  CJ  parameters  and  flow  fields, 
we  expect  that  this  will  be  true,  but  we  have  also  provided  a  further  check  by 
repeating  the  calculations  with  computational  cell  sizes  of  2  x  10  4  m,  1  x  10  4  m 
and  0.5  x  10  4  m.  Table  2  summarises  the  results  of  these  simulations.  As  a 
decrease  in  computational  cell  size  also  directly  decreases  the  time  step,  it  is 
necessary  to  run  for  increasingly  larger  numbers  of  time  steps  as  the  cell  size  is 
reduced  if  the  profiles  and  CJ  parameters  are  to  be  compared  at  the  same  time.  We 
have  used  a  run  with  Ax  =  4  x  10"4  m  after  2000  time  steps  (t  =  105  ps)  for 
comparison  purposes.  As  Table  2  shows,  at  l  =  105  ps  the  detonation  velocity  is 
still  4%  higher  than  the  analytically  calculated  CJ  velocity.  These  results  illustrate 
in  particular  the  large  number  of  time  steps  required  before  the  initial  disturbance 
finally  decays  to  a  steady-state  detonation.  The  CJ  parameters  at  t  =  105  ps  with 
A  x  --  2  x  10'4  m  agree  almost  exactly  with  those  calculated  using  Ax  =  4  x  10'4  m. 
The  differences  are  no  more  than  2%,  except  for  the  von  Neumann  pressure,  where 
the  value  calculated  using  Ax  =  2  x  10-4  m  is  closer  to  the  analytical  value. 


Table  2:  Effect  of  Computational  Cell  Size  H^Oj.Ar  in  ratio  2:1:7 


Ax 

(10-4  m) 

4.0 

2.0 

1.0 

1.0 

0.5 

0.5 

Analytical 

NDUMP 

10 

20 

70 

80 

220 

260 

Number  of 
time  steps 

2  000 

4  000 

8000 

8000 

20000 

20  000 

< 

G 

I 

1600 

1600 

1640 

1660 

1640 

1680 

1542.0 

Pa  (105  Pa) 

12.8 

13.0 

14.4 

14.5 

14.5 

15.8 

12.2 

TCJ(K) 

2350 

2360 

2450 

2500 

2500 

2590 

2298 

pvN  (10s  Pa) 

20.7 

22.0 

22.8 

24.0 

23.2 

22.5 

24.4 

T 

max 

2780 

2780 

2700 

2700 

2500 

2590 

2298 

t  =  105  ns  for  Ax  =  4.0  x  10'4  m,  2.0  x  10"4  m  and  1.0  x  10’4  m,  and  131.3  ns  for 
Ax  =  0.5  x  10  4  m. 

NDUMP  =  3.0  x  106  Joules/cell  in  each  case. 


It  is  of  interest  to  examine  the  profiles  of  the  flow  variables  in  the  case  where 
insufficient  energy  is  deposited  in  the  system  and  the  detonation  dies.  This  is 
illustrated  in  Figure  2,  which  shows  pressure,  temperature,  and  reaction-variable 
profiles  at  two  different  times  for  a  simulation  with  Ax  =  2  x  10  4  m,  DSTRBO  =  3.0 
x  106  J,  and  NDUMP  =  10.  The  pressure  and  temperature  profiles  at  t  =  105  ps 
show  that  the  shock  front  is  located  at  cell  number  600,  while  the  reaction  front  is 
located  near  cell  number  330.  At  t  =  210  ps  the  shock  front  is  located  at  cell 
number  1 100,  while  the  reaction  front  is  located  near  cell  400.  This  separation  of 
the  reaction  front  from  the  shock  front  is  typical  of  an  explosive  event  which  fails  to 
establish  itself  as  a  propagating  detonation.  The  shock  front  velocity  versus  time  is 
also  shown  in  Figure  2. 

Another  interesting  case  is  shown  in  Figure  3.  Parameter  values  used  are 
Ax  =  0.5  x  10  4  m,  DSTRBO  =  3.0  x  106  J,  and  NDUMP  =  260.  The  pressure  and 
reaction  variable  profiles  after  4  000  time  steps  are  typical  of  those  of  a  decaying 
reactive  shock,  and  at  this  stage  a  self-sustaining  detonation  has  not  been  initiated. 
After  a  further  1  200  time  steps  however  the  shock  front  and  reaction  front  have 
become  strongly  coupled  and  an  overdriven  self-sustaining  detonation  has  been 
established.  The  point  at  which  this  occurs  is  evident  in  the  velocity  versus  time 
plot. 


Figure  2;  Pressure,  temperature  and  reaction  variable  profiles  at  t  =  105  (a,  b 

and  e)  and  t  =  210  \is  (c,d  and  f)  for  HjO  lAr.  NDUMP  -  10, 

DSTRBO  =  3.0  x  IQ6  J  and  Ax  =  2.0  x  10 ^  m.  Shock  front  velocity  versus  time  also 
shown  (g). 


Figure  2  ( continued ):  Pressure,  temperature  and  reaction  variable  profiles  at  t  = 
105  (a,  b  and  e)  and  t  =  210  \is(c,d  and  f)  for  H  10  lAr.  NDUMP  =  10, 

DSTRBO  =  3.0  x  106  J  and  Ax  =  2.0  x  10' 4  m.  Shock  front  velocity  versus  time  also 
shown  (g). 
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5.  Results  for  Stoichiometric  #2^2 

The  model  for  H^C^/Ar  just  described  has  recently  been  used  in  a  two-dimensional, 
two-species  reactive  code  to  model  the  transfer  of  detonation  between  gaseous 
explosive  layers  [29-31].  These  simulations  were  performed  to  model  experiments 
conducted  at  the  University  of  Michigan  using  gaseous  mixtures  of  hydrogen  and 
oxygen  of  varying  equivalence  ratios  [32,  33].  While  the  simulations  provided 
good  qualitative  agreement  with  the  experimental  results,  they  were  unable  to 
provide  a  quantitative  comparison  because  of  differences  between  the  explosive 
mixtures  used  in  the  experiments  and  in  the  simulations.  In  order  to  better 
reproduce  the  experimental  results,  we  have  slightly  altered  the  model  for  Yi.JO.JAr 
which  has  just  been  described  so  that  it  can  simulate  detonation  in  a  stoichiometric 
mixture  of  gaseous  hydrogen  and  oxygen.  In  this  section  we  describe  how  the 
constants  for  the  model  are  derived  ana  then  display  typical  results  of  flow-field 
simulations.  We  also  describe  a  slightly  different  method  of  initiation  which 
appears  to  be  more  efficient  than  the  one  used  to  initiate  the  HJOJAr  mixture. 

The  important  quantities  to  be  determined  are  the  rate  of  energy  release  and 
induction  time  as  a  function  of  pressure  and  temperature,  and  the  molecular  weight, 
gamma,  and  total  amount  of  energy  released.  The  induction  time  and  rate  of  energy 
release  were  calculated  using  the  CHEMOD  code  and  the  full  set  of  elementary 
reaction  steps  for  the  hydrogen-oxygen  reaction.  Details  of  the  integration  routines 
used  in  CHEMOD  arc  described  in  reference  [  17],  A  CJ  detonation  in  a 
stoichiometric  mixture  of  hydrogen  and  oxygen  at  an  initial  pressure  of  one 
atmosphere  and  a  temperature  of  300  K  has  a  Mach  number  of  5.2822  (determined 
from  the  Gordon-McBridc  code  [34]).  This  allows  us  to  calculate  a  von  Neumann 
pressure  and  temperature  of  32.4  x  105  Pa  and  1822  K  respectively.  The  CHEMOD 
code  was  therefore  run  over  a  range  of  initial  pressures  and  temperatures  around 
these  values.  Runs  were  made  for  initial  pressures  of  5,  20,  30  and  50  x  10"’  Pa  for 
initial  temperatures  between  1500  K  and  2000  K.  We  found  the  time  dependence  of 
the  energetics  of  the  stoichiometric  mixture  to  be  much  simpler  than  the  argon 
diluted  mixture,  and  virtually  independent  of  initial  temperature  and  pressure.  Each 
of  the  temperature  versus  time  curves  in  the  region  of  interest  was  characterised  by 
a  constant  induction  time  of  approximately  2  x  10  7  s,  followed  by  a  rapid  increase 
in  temperature  until  approximately  90%  of  the  energy  had  been  released,  then  a 
slower  increase  until  levelling  off  at  a  constant  value.  The  lime  interval  between  the 
initial  rapid  rate  of  temperature  rise  and  the  final  constant  temperature  was  found  to 
be  approximately  constant  and  to  have  a  value  of  5  x  10  7  s.  This  lime  dependence 
is  well  represented  by  using  a  reaction  rate  law  of  the  form 


1  dw  1.0  (27) 

W  Hi  5  X  10"7  s 


coupled  to  a  constant  induction  time  of  2  x  10'7s. 
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Figure  3:  Pressure  and  reaction  variable  profiles  att  =  4  000  (a  and  b)  and 
t  =  8  000  (c  and  d)  time  steps  for  HfOfAr.  NDUMP  =  260,  DSTRBO  =  3.0  x  106 
J  and  A x  =  0.5  x  l  O'4  m.  Shock  front  velocity  versus  time  also  shown  (e).  These 
diagrams  illustrate  the  formation  of  an  overdriven  detonation  from  an  initially 
decaying  reactive  shock. 
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Figure  3  (continued):  Pressure  and  reaction  variable  profiles  at  t  =  4  000  (a  and  b) 
and  t  =  8  000  (c  and  d)  time  steps  for  IF  OfAr.  SDUMP  =  2bO,  DSTRBO  =  3.0  x 
10 6 ./  and  At  =  0.5  x  I0'4  rn.  Shock  front  velocity  versus  time  also  shown  (e).  These 
diagrams  illustrate  the  formation  of  an  overdriven  detonation  from  an  initially 
decaying  reactive  shock. 


Three  quantities  remain  to  be  determined;  raw,  y,  and  Af.  Rather  than  try  to 
estimate  these  from  averages  over  appropriate  values  lor  reactants  or  products  as  we 
did  lor  M:/OVAr,  we  adopted  a  different  approach.  Our  primary  objective  in  the 
development  of  this  model  for  stoichiometric  H-,/0,  was  to  accurately  simulate  a 
steadily  propagating  CJ  detonation.  In  particular,  we  required  our  simulation  to 
accurately  reproduce  the  CJ  velocity,  pressure  and  temperature.  Hence  we  treated 
raw,  y,  and  A E  as  variables  (within  appropriate  limits),  and  used  analytical 
expressions  to  find  values  for  these  variables  which  reproduce  the  CJ  parameters. 
For  stoichiometric  H2/02  the  Gordon-McBridc  code  calculates  the  CJ  parameters  to 
be  V  tV  =  2800  m/s,  PCJ  =  18.0  x  105  Pa,  and  TCJ  =  3600  K.  We  were  unable  to 
reproduce  these  values  exactly  by  varying  mw,  y,  and  A E,  but  found  tiiat  the  choice 
of  mw  =  13.2,  y  =  1.37,  A£  =  4.2  x  106J/kggave  analytical  values  of  \'CJ  = 

2800  m/s,  PCJ  =  18.0  x  105  Pa,  and  Ta  =  3200  K. 

In  modelling  the  experimental  results  the  more  important  variables  to  simulate 
accurately  are  the  detonation  velocity  and  pressure,  and  we  expect  the  1 1  % 
discrepancy  in  CJ  temperature  to  have  little  effect  on  the  value  of  these  variables. 
Hence  we  have  adopted  the  parameter  values  described  above,  and  collected  them 
together  in  Table  3  for  reference. 


Following  discussions  with  K.  Kailasanath  of  the  Naval  Research  Laboratory,  we 
have  also  used  a  different  method  of  initiation  for  stoichiometric  H2/02.  The 
procedure  w'e  follow  is  to  first  calculate  the  initial  shock  pressure  and  temperature 
assuming  that  the  material  is  non  reactive;  these  are  the  von  Neumann  values, 
which  we  have  already  calculated  to  be  PVN  -  32.4  x  105  Pa  and  TVN  =  1822  K. 
Rather  than  deposit  a  fairly  large  amount  of  energy  into  NDUMP  cells  at  t  =  0  we 
now  simply  set  the  pressure  and  temperature  to  the  von  Neumann  values.  Even 
with  the  initial  fluid  velocity  in  these  cells  set  to  zero  this  is  sufficient  to  initiate  a 
stable,  propagating  detonation.  The  advantage  of  this  method  is  that  less  excess 
energy  is  deposited  initially  so  the  system  approaches  the  steady  state  more  quickly. 
Another  advantage  is  that  it  removes  the  unphysically  high  temperature  which 
persists  at  the  boundary  long  after  the  initial  deposition  of  energy. 


Table  3:  Parameter  Values  for  Stoichiometric  H 

7 

Induction  time  :  x  =  2  x  10  s 

Reaction  Rate  law  :  1/w  dw/dt  =  -  1.0  /  (5  x  10  7  s) 

mw  =  13.2 

Y  =  L37 

AE  =  4.2  x  106  Joules/kg 

Analytical  CM  \ attics  from  above  constants: 

V  =  2800  m/s 

P  =  18.0  x  105  Pa 

Tu  =  3200  K 


Figure  4  shows  pressure,  temperature  and  velocity  profiles  for  stoichiometric  H2/02 
alter  20  (XX)  time  steps,  or  t  =  560  jas.  We  used  Ax  =  4  x  10  4  m  and  NDUMP  =  3. 
The  shock  velocity  versus  time  is  also  shown.  The  pressure  and  velocity  profiles 
are  similar  to  those  for  H^O^Ar  in  Figure  1  and  show  the  same  self-similarity 
behaviour.  The  temperature  profile  is  also  similar,  but  the  different  method  of 
initiation  has  resulted  in  much  more  realistic  behaviour  at  the  left  boundary  where 
the  initial  disturbance  was  applied.  Also  shown  in  Figure  4  is  the  detonation 
velocity  versus  time,  which  shows  that  the  velocity  has  reached  a  steady  state  value 
of  approximately  2807  m/s  after  only  4000  time  steps.  The  CJ  pressure  estimated 
from  Figure  4  is  18.0  x  105  Pa  and  the  CJ  temperature  is  3. 00  K.  The  analytically 
calculated  values  for  the  parameters  used  in  the  model  are  2800  m/s,  18.0  x  105  Pa 
and  3200  K  respectively.  The  peak  pressure  of  29.2  x  105  Pa  is  10%  lower  than  the 
von  Neumann  value  of  32.4  x  105  Pa,  and  there  is  still  a  large  temperature  overshoot 
of  approximately  10%. 
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Figure  4 :  Pressure,  temperature  and  particle  velocity  profiles  for  stoichiometric 
HfO-,  at  t  =  560  |i.v.  Ax=  4.0  x  10 4  m  and  NDUMP  =  3.  Detonation  velocity 
versus  time  also  shown.  These  diagrams  are  similar  to  those  in  Figure  I  for 
HjOfAr,  but  in  this  case  both  the  induction  time  and  rate  of  energy  release  is 
constant.  A  different  method  of  initiation  has  also  been  used. 


ressure.  xIO  1  iPa) 


Figure  5:  Pressure,  temperature  and  particle  velocity  profiles  for  stoichiometric 
H2i02  at  t  =  70  j-Uv.  Ax  =  0.5  x  10 4  m  and  NDUMP  =  25.  Detonation  velocity 
versus  time  also  shown.  These  diagrams  show  the  same  flow  variables  as  the 
previous  Figure,  but  in  this  case  there  are  eight  times  as  many  computational  cells 
in  the  reaction  zone. 
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We  have  repeated  these  simulations  using  computational  cell  sizes  of  2  x  10'4  m, 
1  x  ICC*  m,  and  0.5  x  10‘4  m  to  again  check  for  any  dependence  of  the  results  on  cell 
size.  Figure  5  shows  pressure,  temperature  and  velocity  profiles  calculated  using  A* 
=  0.5  x  10-4  m  after  20  000  time  steps,  corresponding  to  t  =  70  ps.  The  profiles  are 
very  similar  to  the  ones  calculated  using  Ax  =  4  x  10-4  m,  the  only  difference  occurs 
in  the  peak  values  of  the  variables.  The  peak  pressure  is  now  approximately  3 1 .5  x 
105  Pa,  or  only  3%  lower  than  the  von  Neumann  value,  and  the  temperature 
overshoot  has  been  reduced  to  approximately  2%.  The  detonation  velocity  versus 
time  curve  again  shows  the  velocity  approaching  a  value  of  2807  m/s.  We  used 
NDUMP  =  25  for  this  run  and  the  detonation  velocity  curve  shows  that  this  was  just 
sufficient  to  initiate  detonation.  The  thickness  of  the  curve  is  due  to  a  very  rapid 
oscillation  which  is  too  fine  to  resolve  on  the  scale  shown.  This  oscillation  is  not 
physical  but  is  caused  by  the  "step  like"  position  versus  time  profile  displayed  by 
the  shock  front.  The  Courant  condition  forces  the  shock  to  take  approximately  five 
or  six  time  steps  to  traverse  a  computational  cell,  which  means  that  the  velocity 
continuously  decreases  for  five  or  six  steps  and  then  instantly  jumps  back  to  a  value 
slightly  below  the  value  it  had  at  the  start  of  the  cycle.  These  oscillations  are  also 
present  in  the  detonation  velocity  plot  in  Figure  1,  but  are  too  fine  to  be  seen  on  that 
scale. 

The  reduction  in  the  temperature  overshoot  which  has  been  achieved  by  using  a 
much  smaller  computational  cell  size  can  also  be  achieved  with  larger  cell  sizes  if 
the  reaction  zone  is  limited  to  a  minimum  width.  The  temperature  profile  in  Figure 
6a  has  been  calculated  using  Ax  =  4.0  x  10'4  m  and  a  temperature  correction  of  the 
form 


7j  -  min  ( 1  >  ri<  Ti  ►!  > 


(28) 


As  can  be  seen  by  comparison  with  Figure  4,  the  use  of  equation  (28)  has  reduced 
the  overshoot  from  8.8%  to  3.8%.  The  temperature  profile  in  Figure  6b  is  also  for 
stoichiometric  HJ02  with  Ax  =  4.0  x  10  4  m  but  using  a  stronger  limitation  of  the 
form 


7j  -  min(rI.2,r1.1,7-(1r1,1,r1^) 


(29) 


This  reduces  the  overshoot  to  less  than  1%.  Figure  7  shows  the  result  of  applying 
equation  (28)  to  Hj/O^Ar  with  Ax  =  4.0  x  10'4  m.  Comparison  with  figure  1  again 
shows  a  reduction  in  overshoot  from  around  20%  to  approximately  1%.  It  is 
important  to  note  that  use  of  equation  (28)  on  (29)  does  not  introduce  any  changes 
in  the  other  variables  in  the  simulation,  and  the  only  effect  on  the  temperature 
profile  is  the  significant  reduction  in  the  magnitude  of  the  temperature  overshoot. 

We  have  also  investigated  the  effect  of  the  number  of  cells  used  for  initiation, 
NDUMP,  and  the  total  number  of  time  steps,  MAXSTP,  on  our  estimates  of  the  CJ 
parameters.  The  results  are  shown  in  Table  4  for  the  fixed  computational  cell  size 
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of  Ax  =  4  x  10'4  m,  We  used  the  induction  time  and  reaction  rate  expression 
derived  for  stoichiometric  Hj/Oj,  but  slightly  different  values  for  y,  mw  and  A E. 

The  results  are  as  expected,  increasing  MAXSTP  and  decreasing  NDUMP  lead  to 
increasingly  more  accurate  estimates  of  the  CJ  parameters.  The  table  does  illustrate 
the  asymptotic  nature  of  the  converge,  with  very  large  increases  in  MAXSTP 
required  to  gain  very  small  decreases  in  the  percentage  error. 


' ;  ’  ;  .  •  3 '  z  r  a  •  :  e  - u  r  3  e 1 


CompL.fat.onai  ceC  number 


Figure  6:  Temperature  profiles  for  stoichiometric  H2I02.  Ax  =  4.0  x  10^  m, 
NDUMP  =  3  and  ISTEP  =  20  000.  Reaction  zone  limited  to  a  minimum  width 
using  Ti  =  min  (TiJf  T  ,  Ti+])  in  Figure  6a,  and  T.  =  min  (Ti  r  Ti  r  T,  Ti+J,  Ti+2)  in 
Figure  6b.  The  reduction  in  temperature  overshoot  obtained  by  these  prescriptions 
is  clearly  evident. 
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Figure  7:  Temperature  profile  for  HfOfAr.  Ax  =  4.0  x  /0  4  m,  NDUMP  =  10, 

’  DSTRBO  =  3.0  x  I06  J,  /STEP  =  20  OOO.  Reaction  zone  limited  to  a  minimum 

width  using  r  =  min  (Tj  r  I  ,  T+/).  Comparison  w  ith  the  temperature  profile  in 
Figure  I  shows  that  the  above  simple  limiting  procedure  has  completely  eliminated 
the  temperature  overshot. 


Table  4:  Effect  of  MAXSTP  and  NDUMP  on  accuracy  of  estimates  of  CJ  parameters 


MAXSTP 

1  000 

8000 

20  000 

1  000 

20  000 

Analytical 

NDUMP 

10 

10 

10 

5 

5 

Va  (m/s) 

2  800 

2  680 

2  666 

2  700 

2  650 

2  650 

Pa  (105  Pa) 

23.5 

20.0 

19.0 

20.5 

18.1 

18.4 

TCJ(K) 

3  500 

3  300 

3  280 

3  300 

3  230 

3  245 

NOTE:  Second  method  of  initiation  used 

7=  1.35,  mw=  15.0,  AE  =  4.0x  106 
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6.  Conclusion 


The  results  discussed  in  the  previous  section  demonstrate  the  effectiveness  of  the 
computer  code  described  in  this  report.  The  code  has  been  used  to  simulate 
steadily  propagating  detonations  in  stoichiometric  mixtures  of  H2  and  02,  and  also 
H2/02  mixtures  diluted  with  Ar.  The  results  of  the  simulations  agree  well  with 
variable  profiles  calculated  from  solutions  using  self-similarity  analysis,  and  the  CJ 
parameters  estimated  from  the  profiles  agree  well  with  values  calculated  from 
analytical  expressions.  The  CJ  parameters  are  independent  of  the  computational 
cell  size,  and  both  the  undershoot  in  calculated  peak  pressure  and  overshoot  in 
calculated  peak  temperature  decrease  with  increasing  grid  resolution  in  the  reaction 
zone.  The  effect  of  the  method  of  initiation,  and  number  of  time  steps  included  in 
the  simulation,  on  the  accuracy  of  the  results  has  also  been  investigated. 

We  anticipate  that  replacement  of  the  polytropic  equation  of  state  by  the  HOM 
equation  of  state,  described  by  Mader  [  1  ],  and  replacement  of  the  first  order 
Arrhenius  kinetics  by  a  more  sophisticated  description  of  the  material 
decomposition,  such  as  the  one  described  by  Kim  and  Sohn  [6],  will  enable  particle 
size  effects  in  condensed  phase  heterogeneous  explosives  to  be  investigated. 
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